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ITERATIVE  METHODS  FOR  ELLIPTIC  PROBLEMS  ON  REGIONS  PARTITIONED  INTO 
SUBSTRUCTURES  AND  THE  BIHARMONIC  DIRICHLET  PROBLEM 

Olof  B.  Widlund 

Courant  Institute  of  Mathematical  Sciences 

New  York  University 

251  Mercer  Street,  New  York,  N.Y.  10012 

U.S.A. 


A  finite  element  problem  can  often  naturally  be  partitioned 
into  subproblems  which  correspond  to  subregions  into  which 
the  region  has  been  partitioned  or  from  which  it  was  origi- 
nally assembled.  A  class  of  methods  are  discussed  in  which 
these  subproblems  are  solved  by  direct  methods  while  the 
interaction  between  the  subregions  is  handled  by  a  conjugate 
gradient  algorithm.  The  mathematical  framework  for  this  work 
is  provided  by  regularity  theory  for  elliptic  finite  element 
equations  and  block-Gaussian  elimination.  The  same  tools  can 
be  used  to  provide  a  reinterpretation  of  Bj^rstad's  work  on 
very  fast  methods  for  the  biharmonic  Dirichlet  problem  on 
rectangles  and  a  variant  of  an  algorithm  introduced  by 
Glowinski  and  Pironneau  for  mixed  finite  element  approxi- 
mations of  the  same  elliptic  problem  on  general  regions.  The 
relationship  between  substructuring  and  finite  element- 
capacitance  matrix  techniques  is  also  briefly  explored. 

INTRODUCTION 

In  this  paper,  we  will  discuss  some  special  iterative  methods  for  elliptic  finite 
element  problems  defined  on  regions  regarded  as  unions  of  subregions.  We  are 
interested,  in  particular,  in  the  design  of  algorithms  for  which  the  interaction 
between  the  subproblems,  i.e.  the  discrete  elliptic  problems  on  the  subregions, 
is  computed  with  the  aid  of  a  conjugate  gradient  method,  while  the  subproblems 
themselves  are  solved  by  a  direct  method.  This  approach  differs  from  the  standard 
one  used  in  industry  where  direct  methods  are  employed  throughout. 

The  partition  of  elliptic  problems  into  subproblems  is  a  very  natural  idea  that  is 
widely  used  in  practice  and  which  has  been  discussed  in  the  engineering  litera- 
ture for  at  least  twenty  years;  see  e.g.  Prezemieniecki   [18,19).  Thus  the 
modeling  of  an  entire  mechanical  system  can  profitably  be  organized  by  discreti- 
zing  the  partial  differential  equations  by  finite  elements  on  subregions.  These 
tasks  and  the  factorization  of  the  resulting  stiffness  matrices  for  the  subprob- 
lems can  be  assigned  to  different  engineering  groups  and  computer  systems  or 
processors  with  coordination  required  only  at  the  interfaces  between  the  sub- 
structures to  assure  matching  of  finite  element  triangulations  and  solutions. 
These  ideas  are  particularly  attractive  for  parallel  computing  and  in  the  case 
when  some  of  the  substructures  are   identical  or  previously  analyzed  as  in  the  case 
when  a  simulation  is  repeated  after  the  redesign  or  damage  of  one  or  a  few  of 
the  substructures.  The  final  phase,  in  which  the  interaction  is  accounted  for, 
does  not  lend  itself  equally  well  to  parallel  computing  and  it  also  represents  a 
significant  fraction  of  the  total  computational  work  even  on  a  sequential  computer. 

Three  alternative  methods  are  discussed  in  this  paper  and  numerical  results  for  a 
model  problem  are  presented.  These  algorithms  are  analyzed  in  detail  in  a  forth- 


coming  paper  by  Bjcrstad  and  Widlund  [4],  in  which  a  complete  theory  for  conform- 
ing Lagrangian  finite  element  approximations  of  second  order  elliptic  proDlems  in 
the  plane  are  given.  The  development  of  this  theory  requires  the  use  of  ellip- 
tic regularity  theory  for  problems  on  Lipschitz  regions;  see  Grisvard  [15],  and 
some  apparently  new  finite  element  results'.  We  regard  the  extension  of  these 
results  to  more  general   elliptic  problems  Si  i  siwifficaot  open  problem. 

Earlier  work  on  algori  thms  of  thii  kind  Is  re|irttV-ted  "in  Concus,     Golub  and  O'Leary 
[6],    who  gave  numerical   resMtS  for  PoiS'Son's  probleni  on  T-shaped  regions.    For  an 
analysis  of  this     method  see  Bjiirstad- and  VfidTund    Ml-     T>'0  interesting  new  pre- 
conditioners  for  the  same  special   prob^em.-are-jintrodiiced  in  the  contribution  to 
these  proceedings  by  Golub  and  Mayers    [12].     In  their  numerical   experiments  a 
comparison  is  made  bet"2en  these  algorithms  jind  the  one  characterized  as  second 
best  in  our  previous  paper;  cf.   Bjiirstad  and  Widlund    (3).     Methods   inspired  by 
Schwarz'     alternating  method  and  by  control   theory  are  considered  in  Glowinski, 
Periaux     and  Dihn   [10].     For  a  discussion  uf  one  of  these,  see  Bjiirstad  and 
Widlund   [4].   There  has  also  been  extensive  activity  in  the  Soviet  Union  in  this 
area.  We  are  not  yet  well  acquainted  with  this  work  and  a:  e  grateful   to  Yu. 
Kuznetsov  and  Maximili.-.n  Dryja  f  or^  drav«  ng^  our  attention  to  it.   After  our  first 
series  of  numerical   experiments,  we  reCfii^ed^a- preprint  of  one  of  Dryja's   papers 
[8],  from  which  we  first  learned  of  what  we -currently  regard  .is   the  best  algorithm. 
In  that  paper  Dryja  analyzes  the  case  of  Laplace's  equation  on  L-shaped  regions. 

The  techniques  devel'-ped  for:  thesubstructuring  problems',  which  combine  block 
Gaussian  elimination  and  regularity  theorems  for  elliptic  problems,  are  also  use- 
ful   in  examining  certain  iterative  methods  for  the  biharmonic  problem.    Even   on 
rectangular  regions,   thi'^,  problem  is  of  considerable  diffi'culty,   since  separation 
of  variables  cannot  be  used  directly.     By  introducing  an -auxiliary  variable,   known 
as  vorticity  in   fluid  dynamics,   the  problem  can  be  reformulated  as  a   second  order 
elliptic  system.     Many  iterative  methods  have  been  developed  to  solve  this  system 
in  which  a  different,   simpler  boundary  value  problem  is  solved  in  each  step;   see 
Bjjirstad    [1,2]     and  Glowinski  and-Pironneau   [11]     and   the  many  references  given  in 
those  papers.      In  this  paper,  w«  introduce  an  alternative  to  an  optimal   method 
giVen  by  Glowinski   and  Pironneau  l^^]     and  give  a  reinterpretation  of  Bj^rstad's 
work. 


While  for  substructured  problems'  we  seek  to""  solve  a   large  linear  system  by  an 
iterative  method  which   involves  smaller  problems,   there  are  occasions  when  we  are 
willing  to  solve  a   larger  linear  system  repeatedly,   to  obtain  the  solution  of  a 
smaller  system.     Such  algorithms  are  known  as  capacitance  matrix  methods;   see 
O'Leary  and  Widlund    [17],   Proskurowski   and  Widlund    [20,21)      for  a  discussion  and 
references  to  the  literature.      In  our  last  section,  we  explore  the  relationship 
between  substructuring  and  finite  element-capacitance  matrix  methods. 

SUBSTRUCTURED  FINITE  ELEMENT  PROBLEMS  AND  BLOCK  FORM  OF  THE  STIFFNESS  MATRICES 

To  simplify  the  discnssion,  we  confine  ourselves  to  problems  defined  on  the  union 
R  of  n,  ,  Oo  and  r,  and   to  the  Dirichlet  case.     Here  n,   and  n-  are  plane,   bounded, 

nonintersecting     regions  and  r,  the  intersection  of  their  closures.       The  bound- 
aries of  Q-,   and  Q-  ^'"^  '"i  '-'  ''3  ^"'l  ""p  '^  ^3  '     '"espectively ,  and  the  boundary  of 
n  is  r,  u  r-.     We  assume   that  these  subregions  are  curvilinear  polygons,   i.e.    in 

particular  they  are  Lipschitz  regions;   see  Grisvard   [13].       A  linear,   second  order, 
positive  definite,  selfadjoint  elliptic  operator  is  defined  on  Q.      Its  syriw^tric 
bilinear  form  is  denoted  by     ao(u,v);   see  Ciarlet   [5]    or  Strang  and  Fix    [22]. 
A  simple  example   is   given  by   tKe  Laplace  operator  with  homogeneous  Dirichlet 
condition   for  wnich 


a„(u.v)  = 


VU-7V    dx 


u,v  e  H-(n) 


Here  HqCq)   is  the  subspace  of- elements,  with  zero  boundary  values,   of  the  Sobolev 
space  of  square  integrahle  fufKtions  with  square  integrable  first  derivatives. 
Triangulations  .of  n-j.  and  Q^  are  intrj>duced  in  such  a  way  that   the  nodes  on  r3 
coincide  and  rj^j.^fotl^s-   .elepienit, boundaries.     We  assume  that  each  degree  of 
freedom  of  the  finite  e'lement  subspa'ces  is  associated  with  a  node  and  a  basis 
function  in  theffin\te,de-l^?n>ent  space  jnd  that- t{ie,. support  of  any  basis  function 
coincides  with  the  iria'ngl.e's  to  whicti  its  node. belongs.   An  element  of  the  stiff- 
ness matrix  hag  the  torm    .sJ^.^ri^)  ,c,)tilneni'.'^.c.ar\(i  <^.  are  basis  functions,  and  it 

therefore  vani'sheS'^'iinles's  'th^  twO  nodes  bel&n|"to  a  common  triangle. 


It  is  easy  to-Siefr^m  the  tfefinitfcst'of  the  bilinear  form    that 


\-^(u.v) 


f.a 


n; 


(u.v) 


(u,v) 


(2.i; 


and  that  therei^ore  the  stiffness  Watrix'of  a  problem  on     n     can  be  constructed 
from  those  of -f!i   and  nn.     The-  same-relat^bn  holds  for  any  pair  of  nonoverlapping 
subregions.     This  fact^  freqi)eht3'y  used  in  practice  to  construct  stiffness 
matrices   from"tb« 'st+ffness  aatHc'es  of*th6?lTidividual    triangles;   see  Strang  and 
Fix   [22J.  ..r  :^-.      •■'•:?'-.  --     -...J  :i  -    :r  ." 


We  write  tbe-'atifefJriessinBtnces:  correspen^jing- tb  Qy  and  n^  respectively  as 
;-3  j'*:r'  ?  'o'-  ^     HC=f*    ,:  -«.'„;.;     ': 


[i; 


;j5Tiit     r/^^fjp--!'    ,;  and. re  K&? 


23 


^^23 

.(2) 

^33 


(2.2) 


:o'. 


where     K, ,  ,    represents  couplings  between, the  pairs  of  nodes  in  n,  ,     K,,  couplings 
between  the  pairs  be^To'nging  to  n.j«'^n'd  f,  respectively  and  K  ,3     couplings  between 


nodes  on  r,,  etc 


Written   in   variational    form  a  Dipl^hlet  problem, on  n,    has  the  form 


vu, -vv,    dx  = 
h       h 


f  V,    dx 

n 


VV^€S 


n  HQ(n^)    , 


u^,  e  s"(n^)    , 


Vh  =  9c 


where  S  (n^ )  c  h  (n)  is  a  finite  element  subspace  and  YqU^  is  the  trace,  i.e. 

the  restriction  of  u.  to  the  boundary.   In  matrix  form  this  problem  can  be   writ- 
ten as 


*"  *^11   '^13 


(2.3) 


where  x-j  and  x,  are  the  vectors  of  nodal  values  corresponding  to  the  open  set  n-, 
2  respectively.   The  components  of  b,  are  constructed  from  the  right  hand 


and 

side  ■'  f  and  the  Dirichlet  values  on  r 


and  b,  is  the  vector  of  Dirichlet 


values  on  r,. 

If  on  the  other  hand  Neumann  data  g^^   are  given  on  r,  ,  while  we  still  have 
Dirichlet  data  g„  on  r, ,  the  problem  is  of -the  form,  -  : 

7u,.7v^  dx  =  I  f  v^  d..^J^v^^ds;i  [y^y^f  "  njln^.r,)  . 


u^,e  S 


'1 


^0  "h  'W 


Here  yq  u^  is  the  trace  on  r^  'and  HQ(n^,r^)-  the  subspace  of  H  (n^ )  with  vanish- 
ing trace  on  r,.  In  this  case  x,  is  a  vec,tor-of  unknpwns  and  the  linear  system 
is 


^11 

^13"      ^^33 


<i3 1  ^^ll..  ["ii 


kUM 


lb  J 


(2.4) 


We  note  that  the  vector    b,- vanislies  ;i?f -f  and  gjj-are  zero, in  which  case  the  vector 
b,  represents   the  Neumanri  d'a.ta,  ;o"h,  Z,^- .  •        '.J— '-u'       jri'   '•■   - 

.-■-1-  '  ""•  '*''^.i 

The  stiffness  matrix  of  thVfihiire  prob3etit-i's  "of  tHeifdrrt, ' 


K  *. 


^n 


0  ,.      .Kj 
.     ^2 


'13       "23 

where,   by  (2.1)   and  a  simple  fornjuit^tion,,,  ^,,^ 


'^23 
^'33 


'^33   "^33   ,33  . 


(2.5) 


We  note  that  the  degrees  of  freedom  have  been  partitioned  into  three  sets  of  which 
the  third,  the  separator  set,  corresponds  to  the  nodes  of  r3.  From  the  point  of 
view  of  graph  theory,  the  undirected  graph  of  K  becomes  disconnected  into  two 
components  if  the  nodes  of  the  separator  set  and  their  incident  edges  are  removed. 
If  conforming  finite  elements  are  used,  then  it  also  follows  from  the  assumptions 
on  a-(u,v)  that  K  is  positive  definite,  symmetric  and  as  a  consequence  so  are 

1 1  *   7? '   ^^     33 ' 
We  consider  the  linear  system  of  algebraic  equations  of  the  form 


Kx  = 


^11 

0 

^13 

''l 

0 

^22 

K23 

"2 

J 

J 

^13 

hi 

^33 

""2 

'^ 


(2.6) 


By  using  bloclc-Gaussian  elimination,  we  can  reduce  this  system  to  the  positive 
definite,   syrmetric  system. 


SX3  = 


^33  ■  *'~13   h]    ^3 


h3  '^22  '^23'''3 


'  ^3  '  ^13  h]   ^]   '   h3  h2  ^2  '  ^3 


(2.7) 


It   is  common  practice  to  complete  the  process  by   solving   (2.7)   by  a  direct  method. 

The  right  hand  side     bj  can  be  obtained  at  the  expense  of  solving  the  two  suborob- 
lems  on  ^1   and inj,  multiplying  the  resulting  vectors  by  the  sparse  matrices   KJ., 
and  K^3         respectively     and  by  subtracting  the  resulting  vectors  from  b,. 
From  now  on,  we  will   consider  only  the  case  when  b-j   and  bp  are  zero.     SuCh  a 
reduction  can  of  course  easily  be  accomplished. 

The  matrix  S,  is  a  so-called  Schur  complement;   see  Cottle   (7),  and  can  be  expen- 
sive to  compute  andr^tore.     However,  we  notice  that  Sy  can  be  computed  for  a  given 
vector  y  at  the  expense  of  solving  the  two  subproblems  with  the  sparse  right  hand 
sides   K,,y  and   K-AV,"  respectively,  and  certain^  sparse  matrix  and   vector  operations. 
In  the  next  section,  we  will   develop  iterative  methods  which  only  require  S   in 
terms  of  such  matrix-vector  products. 


The  cost  of  computing  Sy 'depends  primarily  on  the  efficiency  of  the  solvers  for 
the  subproblems.  It  should  aVso  be  rwted  ttiat  if  a  Gaussian  elimination  method 
is  used,  advantage  can  be  taken  of  the  sparsity  of  the  vectors   K,  ,y  and   Kp-^y. 

Thus  when  the   lower 'tistafignlar-systWiiS-9)f  ec(Uations  are  solved,   the  computation 
can  begin  with  the  first  equation  which,  h^^  |L,.non2erfl  right  hand   side.    Similarly, 
the  solution  of  the  upper  triangular  systems' can  be" stopped  as   soon  as, all    the 
components  of   '^ij,'^|3J^     ^P,d   Kjl  K-^.X^   n.ece^jary  for  computing   kLujI    K-^y)  and 

Kp-jlK^p  '^23^'*  ^*^®  ''®^"  found.     This  can  effectively  reduce     the  size  of  the 
triangular  systems  necessary  to  carry  out  the  iteration  steps.      It  is   thus  parti- 
cularly advantageous   if  all   the  variables  _at  nodes   adjacent  to  r3  are  ordered 
late.      It  should  be  not^d,  however,   that  such  a  constraint  may  be  hard  to   impose 
on  existing  software  or  that  it  may  lead  to  an  increase   in  the  time  required   to 
factor  K,,   and  K.-   into  their  .triangular  factors. 


We  will   need   the  Schur  complements  wi't^  respect  to  the  matrices 
defined  in  (2.2).   They  are> 


.(1) 


and  K 


(2) 


■(1)   -   .(1) 
'         "   ^33 


Using  (2.7)  and 


f  ■  ■■  -1    " 
^13  "^11    '^13 


;2.8),  we  find  that 

S  = 
(1) 


arid 


:(2) 


,(2) 
^33 


1^23  ^22 


^23 


;(!)   ,s(2) 


(2.8) 


(2.9) 
(1). 


The  mappings  S  and  S*  ' '   play  an   important  role   in  what   follows.      The  vector  S*    'y 
can  be  computed  by  solving  a     Dirichlet  problem  (2.3)   and  then  applying  the  matrix 
of   (2.4),  which  corresponds   to  a   Neumann  case,   to  the  solution  vector.      This  can 
be  seen  by  a  straightforward  computation: 


11 
T 

13 


13 

(1) 
33 


,-1 


Ml 


M3 


I 


S'l'y 


From  our  previous  discussion     we  see  that  this   is  a   Dirichlet-Neumann  map  which 
takes  the  Dirichlet  data  y  on     r3   into  the  discrete  Neumann  data  S^^^y    on   the 
same  set.     For  a  discussion  of  the  continuous  case,   see  next  section. 


CONJUGATE  GRADIENT  ALGORITHMS   FOR  SUBSTRUCTURED  PROBLEMS  AND  AN   INFORMAL   THEORY. 

The  general    theory  of  conjugate  methods   is  quite  well    known  and   it  will    therefore 
be  discussed  only  very  briefly;   see  Concus,  Golub  and  O'Leary   [5],   Hestenes    [14] 
or  Luenberger    [16] . 

Let  Ax  =  b     be  a   linear  system  of  algebraic  equations  with  a  posi ti ye  defini te, 
symmetric  matrix  A.     Let  x'°'  be  arr  fnitfar  gueis  4nd  fW  =  b  -  Axl"'     the 
initial    residual.     The  k-th  iterate  in  the  standard  conjugate  gradient  method, 
x^^i,  can  then  be  characterized  as.  the  minimizing  element  for  the  variational 
problem 

''•'miri(^)y'^Ay  -  y'^b'  ■'^' 

where  y-x'°^  varies   in  the  linear  space  spanned  by  r^   ; ,  Ar^    ' A  "  r*    '. 

By  expanding  in  eigenvectors  of  A,  4t'£an  be  established  that 


X 


(M_,\T  nrv(k).,^/,.(0). ;,T  .,   (0) 


is  bounded  from  above  by. 


A(x^'''-x)/{x^"'-x)'   A(x^"'-x) 


,niin:    ,,  max, =:(%,- xp(x))^, -•;.  (3.1) 


■p6P^_.,  xeo(A) 

See  Luenberger   [16].     He'-e.x.is  the  exact- solution:,^  P^.-|   the  space  of  all   poly- 
nomial    of  degree  k-1  £.1:0  o(A)  the's0ectrtim-of  A.     This  bound  can  be  used  to 
establish  that  the  convergence  fs 'rapid   If.  A  is  we'll 'fond^tloned  and   that  the 
rate  of  convergence  can  Ij  bdund^d  uniformly  for  ehtire  famil ies  of  operators   if 
all   the  eigenvalues  fall   in  a  fixed  interval.  -^jriL 

Preconditioned  conjugate  gradtent  methbds  have  been  studiea  extensively  in  recent 
years.     The  idea  goes  back  to  the  mid-fifties;  see  Hestenes   [14].     Let  Aq  be 
another  positive  definite,   synmetric  operator  for  which  it  is   feasible  to  solve 
auxiliary  systems  of  the  form  AQy  =  c  repeatedly  for  different  right  hand  sides. 
In  one  of  the  versions  of  the  method, 'trie  oriiginal   problem  Ax  =  b  is  transformed 

into  AAZ  y  =  b.     The  iterate  y     '   is  sought  as  the  sum  of  an  initial   guess  y^ 
and  a   linear  combination  of  r^^' ,   AA^  r^^' ,. . .  ,(AAq   )   ""  r°.  If         an  appropri- 

ate inner  product  is  used,  a  convenient  recursion  formula  results,   see  e.g. 
Proskurowski   and  Widlund    [21].      Each  step  of  this  algorithm  requires   the  solution 
of  an  auxiliary  linear  system.      It  is   important  to  note  that  the  estimate   (3.1) 
still   holds,   but  that  now  the  eigenvalues  of  AAZ   ,   i.e.     those  of  the  symmetric 
generalized  eigenvalue  problem  A$  =   xAq*,     are  of  relevance  rather  than  those  of 
A.    It  is  also  worth  noting  that  the  estimate   (3.1)   can  be  used  to  show  particular- 
ly rapid  convergence  if  the  eigenvalues  are  clustered. 


For  the  problem  at  hand,  we  first  consider  the  solution  of  equation   (2.7)  without 
preconditioning.     From  (2.9),  we  see  that  Sy  =  S^''y  +  S'^'y.      If,]^  therefore 
at  least  plausible  that  S     will   be  ill   conditioned  if  s(l)   and  S(''  are.   As  shown 
in  section  2,  s'^)   represents  a  Di richlet-Neumann  map  and  therefore   involves  a 
loss  of  a  derivative   in  L2(ro).     Such  a  map  will    have  a   spectral    condition  number 
proportional    to  the  number  of  nodes  on  r3.      In  order  to  clarify   this  point,  we 
consider  the  continuous  case,   leaving  the  details  on  the  finite  element  case  to 
the  forthcoming  paper  by  Bjjirstad     and  Widlund   [4]. 

Thus  consider  two  harmonic  functions  Ui   and  uj  defined  on  n-^   and  Q2  respectively. 

These   functions  vanish  on  r,    and   r^     and   have  the  same   trace  on   r^.    They  can 

1 
therefore  be  combined   to  form     u  e  H-(.i.).     This   function  satisfies 


a^tu.v, 


=    f(v)     ,       V    V    €    HqIp.)     , 


(3.2) 


Sj^lU.V, 


vu,  -wdx  +   I   VU2- 


where  the  linear  functional    f  has  its   support  on  r^.    It   is  easy  to  show  that 

••vvdx  =  |     l^vds-j     |^vds  =  j!|^]vds. 

"l  "Z  .  r  •    ^Xr  '"S  '"S 

where  v  is  the  nqrjn^][^Citvi^^^  ^(pi  rj^'p^t'.lo  n^.     We  can  therefore  rewrite  equa- 
tion (3.2)  as   .,.""'■' 


-I 


•6. 


■.Vi^r  i, 


\^\''-'f    . 


'3v 


„here   [— ]     corresponds  to  S^.     Following.  l,ions  and  Magenes    [15]   we  see  that  for 
u  e  Hj(n).  ■  ',   ,     .. 


Vq^^^^TO) 


[u  ^mY^(tI)   ; 


Ip  ui 


L^(r3) 


»} 


where   o     is   the  distance  of  a  point  on  r,  to   its  end   points.      It  then   follows, 

from  a  standard  variational   argument,   that     -^ ,    JT  ^'^'^    W      belong^to  the 

1  /  2      *" 
dual   space  of  Hqq 


■^  '"^(r,)';'  i.*-'.  "a  derivative  is, Tost  in  comparison  with  -.q  u. 


'  -  )r.  .5  , 


In  view  of  what  we  have  just  learned,   it  is  natural    to  try  to  find  a   precondi- 
tioner  which.'al^so  invol  veT- tWe'^ess'flf  a^flerfvati^i  ■  in  L2(r3).   A  natural   choice 
would  be  a  tanMntJa-]  .d€rivataye-but_that  is  not  ^a  ^syrrmetric  operator.    Instead  we 
can  use  the  sauafi^X|i?x)t  p/^  ttie"  negative 'of , a  disi^re^tization  of  the  Laplacian  on 
r3.     Such  a  miettipxt  If  practifal   fo*  problei^s.-.tp.  the  plane  and  has  been  tested; 
see  section  4.   We  denote  this  opgrator  as^J.v 

An  even  better  m.eth^,^  ipvg,]veSjJ;he  ^9lution|.Qf  a  system, 
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can  be  accomplished     by  solving  equation   (2.4 
d  the 
pping 

0 
K 


This  solution 

side   (0,y)'     and  then  the  discrete  Dirichlet  problem  on  rj2 

The  relevant  mapping     is  now     SS^)'       since 


M3 
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with  the  right  hand 
using  X3  as  data. 


SS 


:i) 


-1 


It  can  be  shown  that  this  operator  is  uniformly  well  conditioned;  see  Bjjirstad  and 
Widlund  [4] . 

We  note  that  we  need  not  construct  any  auxiliary  operator  when  this  precondi- 
tioner  is  used.   It  is  also  interesting  to  note  that  if  a  rectangular  region  is 
cut  in  half,  and  treated  fully  symmetrical ly,  then  S  =  2S^''  and  the  conjugate 
gradient  iteration  converges  in  one  step. 


NUMERICAL  EXPERIMENTS  WITH  A  SUBSTRUCTURED  PROBLEM. 

We  will  only  gi*e  a  brief  report  on  a  few  of  the  many  experiments  which  we  have 
carried  out  for  the  five  point  approximation  of  Poisson's  equation  on  regions 
which  are  unions  of  two  rectangles.  The. choice  of  such  regions  greatly  simplifies 
the  experiments  and  makes  it  feasible  to  conduct  n.^^ny  experiments  with  very  many 
degrees  of  freedom  since  fast  Poisson  solvers  can  be  used  to  solve  the  subprob- 
lems.  We  note  that  this  simple  finite  difference  approximation  can  be  viewed  as 
a  conforming  finite  element  approximation  oo  a. mesh  of  right  triangles  using 
piecewise  linear  functions;  see  Strang  and  Fix  [22]. 

In  the  experiments  considered  here,  we  consider  the  union  of  two  rectangles  with 
corners  at  the  points  (0,0).  (1.0),  (hl/2),  (0.1/2)  and  (1/8,1/2),  (5/8,1/2), 
(5/8,1),  (1/8,1)  respectively.  The  relevant  parameters  are  the  number  of  mesh 
points  q  on  r3  ,  the  interval  between  (1/8,1/2)  and  (5/8,1/2),  and  the  total 
number  of  degrees  of  freedom  H^   have  used  data  which'  are  consistent  with  an 
exact  solution  u(x,y)  =  x^+y'  .  xe^cos  y.  We  have  found  no  real  difference 
between  the  performance  of  our  method  for  this  and  other  cases. 


We  first  show,  in  Table  I,  how  the  number  .of  iterations  grows  as  a  function  of  q 
using  the  operator  SSd)-'.  We,  stopped. the  iterations  at  the  level  of  the  trunca- 
tion error.  The  initial  gues*  was  the  zero  function.  Ue  note  that  the  overall 
number  of  degrees  of  freedom  increases  quadratical 1y  with  q  and  equals  48641  for 
q  =  127.  .'-■'  r'iwoi:o-  s"-  r^c  '.  ■  •  r.  ■<■•   ,  r- •  : ' 

■    ■  .r  ('  ,5? 


.        .-.      :TABLE,r      '- 

'.'•'■■  -■'^' 

q 

Number  of 
-     iteratiorrs 

Maximum  Error 
"  ■    ■            on  n 

3 
7 

.      -    =  -   -  2    - 

3.66x10"'* 
9.59x10"^ 

15 

.:   n'     .    =  ..-.-  3 

2.45x10'^ 

31 

63 

•   -    ..„ ,  4 
■   -4 

6.09x10'^ 
1.49x10'^ 

127 

5 

3. 02x10°^ 

In  Table  II,  we  compare  this  precondi tioner,  S     with  J,  which  is  the  square 
root  of  the  negative  of  the  discrete  Laplactan  on  r3  ,  for  the  case  when  q  =  127. 


TABLE   II 

No.   of   Iterations 

Precondi tioner:   S 

Precondi tioner:     J 

0 
1 
2 
3 

4 
5 
6 

7 
8 

3.79x10-!! 
1.25x10-2 
7.48x10-^ 
2.56.10-5 
4.42x10"' 
3.02x10-; 
3.02x10-; 
3. 02-10-' 
3.02x10-' 

3.79x10-^ 
3.22x10-2 
4.01xl0--' 
5.26.10-5 
8.74.10-= 
1.05x10-5 
1.33.10-° 
3.08x10-' 
3.03x10-' 

Ten  to  twelve  iterations  were  required  to  decrease  the  maximum  error  by  a  factor 
of  10  if  no  precondi tioner  was  used. 

(1)-^ 
We  also  give  some  information  on  the  spectrum  of  SS     for  q  =  53.  The  small- 
est eigenvalues  were  found,  to.  be  >i  =  1.713,  /j  "  1-777,  .  . . ,  x^  =  1 .992  with 
all  eigenvalues  less  than  or  equal  to  2.000.   The  corresponding  spectrum  for 
SJ"^  shows  a  somewhat  IgS'S  pronounced  cluster  with  x,  =  1.733,  X-,  =  1.804,  ..., 
Xc  =  2.008 '^,','-''2.'628.   '  ^ 

DIRICHLET'S  PROBL^  FOR"YhI  8'IHARMONIC  EQUATION. 

*■   :  ■      ■:  ti- 

The  bihamionic  Dirichlet  problem  has  the  form, 

" '     '      ''!:■"        ^-*  '/     ''"  "•- ■ 

where     Yq'jj  and     y^Ai     are  traces;  t  .e;"the  r9'S:trictions  of    i    and   at/Bv     to   the 
boundary.      For  a     discussion  of  the  existence  of  traces,   a  definition  of  the 
Sobolev  spaces  used;in  this  section  and  a  gensral    introduction  to  elliptic  prob- 
lems;  see  Lions;  and' MageneS   [>5J.'  We-  assume  that  the  region    n    is  plane,  bounded 
and  simply  connected   ;  and"'denotetth€'b(airKlffry -by  r. 

In  our  discussion,  we  will   use  the  following  regularity  result;   cf.   Glowinski 
and  Pironneau    [11]  : 

Let  r     be  sufficiently  smooth,  arid  "assume  that  f  s  L  (n),  gg  e  H  '    (r)  and     ,,. 
g.|  e  H''^(r).-- tterv-t+te-so+trt+en-of -f5^H-s«-ti^'fts3  li- e  w'^X'^)     and  >g(ai06H'  'v). 

It   is  easy  to  show  that  we  owly  .rwed-'to  consider  the  case  where  f  =   0  and  gn  =   0 
in  what  follows..-.-fo-Uowing-&k)w^nsV4-dfKt -Pironneau    [11),  we  attempt   to  solve 
equation   (5.1)   by  finding   x,   th^  trace  of  the  vorticity  <^  =  -a*.   for  which 
Y.|*  =  g.| .     The  relation  between  .Yqu  and  y^*     is  given  by  the  system, 

;  -fi(u  =0     in  n, 

■■•^0"=   ^  '  (5.2) 

-A0  =01     inn, 

Yq*  =   0 

-1 17. 
which  can  be  shown  to  define  a  selfadjoint,  H    (r)-elliptic  operator  A  which 

maps  H"''2(r)  onto  H'''(r).  At  the  expense  of  solving  two  Dirichlet  problems  for 

the  simplest  second  order  elliptic  problem,  we  can  thus  compute  Ax  for  a  given 

X  and  can  then  use  the  conjugate  gradient  method  to  find  x  and  the  solution  of 

problem  (5.1 ). 

There  have  been  many  attempts  to  use  related  methods  for  finite  difference  approx- 
imations of  (5.1);  see  Bjbrstad  (1,2)  and  Glowinski  and  Pironneau  [11]  for  a 
discussion  and  references.  The  idea  is  systematically  developed  by  Glowinski  and 
Pironneau  in  the  framework  of  mixed  finite  elements.  These  have  proven  to  be 
the  most  successful  class  of  methods  for  solving  problem  (5.1).  We  review  some 
of  their  work  stressing  the  similarity  with  the  work  discussed  in  sections  2 
and  3. 

The  mixed  finite  element  methods  are  defined  by, 


Vii). -vv,  dx  =  0  , 


^  %  ^  ^0 


and 


vu,.  .7v.  dx 
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g,u^,  ds  .  vu^  e  M^ 


Here  S  c  H  (n)  is  the  finite  element  subspace, 


^0. 


space  spanned  by  the  basis  functions  associated  with  the  modes 


Hjla) 


and  Mu 


the 

on  r.     f^h  "  ^'^ 
thus  a  subspace  of  S"  which  is  complem«ntary  tcSg.     With  ^^^   denoting       the  vector 
of  nodal   values  of  j^,   in     n,  u-,    the  vector  of  the  negative  of  the  nodal   values  of 
uu  in     n     and     uj   the  vector  values  o^-i^  O"   r,  we  obtain  a   linear  system  of 
equations  with  a  symnetric,  -indefinite  coefficient  matrix, 


0 

K,f      •K,2^ 

^♦1 

0 

•^11 

"1 

X 

b 

(5.3) 


Here  K, ,    is  the  stiffness  submatrix  which  represents   the  couplings  between  pairs 
of  nodfes  in  n     while     Kjg  ""epresents  couplings  between       n  and  r.    The  matrices 
M.  .      .    •  ...      ,         . 


11 


M,  2  and  M-o  are  ma 


ss  submatrices  with  elements  J„  t^'t^   dx  .  The  matrix 


and  the  mass  matrix  ire   positive  definite,  synmetric. 


Diagonal  subsystems 


associated  with  the  vectors  ii  and 


can 


5  ^^■,  diiu  u^   V.OM  be  solved  at  the  expense  of  two 
discrete  problems  for  the  Laplacian.  As  in  the  continuous  case,  we  can  reduce 
the  problem  (5.3)  to  a  linear  system  for  the  trace  of  the  vorticity.  Using  block 
Gaussian  elimination,  we  can- regard  -the -resul  ting  matrix,  A^  ,  as  a  Schur  comple- 
ment. The  syimietry  of  Af,  follows  by  i/i spec t ion.  The  fact  that  it  is  positive 
definite  can  be  shown  by  using  j  result -given  in  Cott'c  [7].  He  uses  Sylvester's 
theorem  to  relate  the  inertia  of  a  block  matrix  with  that  of  a  principal  minor 
and  the  corresponding  Schur  complement.- 


The  matrix  A^  is  an  ap 
used  techniques  from  m, 
region,  the  condition 
This  is  not  surprising 
is  gained.  They  theref 
a  discretization  of  th 
in  the  sense  that  the 
degrees  of  freedom  of 
cally  the  same  as  the 


oroximation  of  the 
athematical  analys 
number  of  A^  grows 
si  nee  A  maps  H"  '' 
ore  suggested  a  pr^ 
e  operator  -(d^/ds 
rate  of  convergenc 
the  discrete  mode" 
inverse  of  the  ope 


operator  A.  Glowinski  and  Pironneau 
is  to  prove  that  in  the  case  of  a  convex 
linearly  with  the  number  of  nodes  on  r. 
2(r)  into  H'''(r),  i.e.  a  derivative 

conditioner,  based  on  the  square  root  of 
"^l*  1.  This  provides  an  optimal  method 
e  becomes  independent  of  the  number  of 

This  precondi tioner  is  of  course  basi- 
rator  J  described  in  section  3. 


In  view  of  the  success  of  the  precondi tioner  S 
we  have  initiated  numerical  experiments 
as  a  preconditioner  for  A^.  Given  a  vec 
the  same  discrete  Laplace  problem  which 
Since  the  continuous  as  well  as  the  disc 
one  dimensional  null  space  of  constants 


(1) 


for  the  substructured  problems, 


which  will  use  the  Di richlet-Neumann  map 

tor  y,  s'''y  can  be  computed  by  solving 

is  required   for  the  computation  of  A.y. 


modified  by  adding  cee    ,  where  e' 
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in  this  case,  the  operator  S 
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Fast  algorithms  for  the  biharmonic  problem  on  rectangles  were  considered  in  the 
1980  Stanford  dissertation  of  Peter  Bjeirstad  [1];  see  also  Bjjirstad  (21.  Although 
the  region  is  quite  special  such  algorithms  are  nevertheless  very  useful  in  a 
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number  of  applications.  The  discrete  model  is  the  standard  13-point  difference 
approximation  of  problem  (5.1)  on  uniform  meshes,  but  to  simplify  the  notations 
we  will    instead  describe  Bjibrstad's  algorithm  for     the  continuous  case. 

The  task  is   thus   the  solution  of  equation   (5.1)  with  a   rectangular  region  ::.    The 
boundary  r  =   Fh  u  fw    ,  where  r^   is  the  union  of  the  two  horizontal   sides.   Because 

of  the  simple  geometry  there  are  alternatives   to  equation   (5.2)  when  we  search 
for  problems  wtiich  ,are  eajy  to  solve.    -Thus  we  can  solve. 


(5.4) 


by  separation  of  the  variables;   ''W^ido 'tMs  by  expanding  in  Fourier  series   in  the 
horizontal   directiDn:and  ■sorvJng'the'>r^J11;i.ng  boundary  value  problems  for  a   set 
of  fourth     order  ordinary  differential  -equations.     We  note  that  the  data   for  equa- 
tion  (5.4)   is  quite' special ,   bCit  that  we  can  always  reduce  our  problem  to  this 
form  by  solving  a  preliminary  separable  problem.      In  the  discrete  case  the  appro- 
priate algorithmic' tools  are  the  fafst  Foi/rier  transform  and  the  Choleski       algo- 
rithm for  five  diagonal    linear  systems  of  equations. 

Equation   (5.4)   can  be  used  to  define  a  mapping-  from  yqw  to     y]*  on  r,.     To  find  a 
preconditioner,  wer.'consider  another*  sepafabl^ 'problem  for  which  we  can  use  a 
Fourier  expansion  in  the  vertical 'di  rection.  ' 


(5.5) 


By  evaluating  y„w  on  T„-for  th#"sdfLitio"n  "of  ,(5.5-) ,' we  define  the  precondi  tioner 
which  has  proven   to  be. highly     Successful >■'  in  ejiirstad's  work.     A     tight  uniform 
bound  on  the  spectrum  of  the  precdriditi dried  operator   is  rigorously  established 
in  his  work.      In  computational   practice,   the  error  decreases  by  a   factor  of  about 
10  in  each  iteration  step.     We  also  note  that  advantage  can  be  taken  of  the 
sparsity  of  data  of  the  problem  and  the  fact  that  in  the  iterations,   the  solution 
is  only  required  on  and  close  to  T\j. 

FINITE  ELEMENT-CAPACITANCE  MATRIX  METHODS 

There  are  occasions  when  a   linear  system  of  equations  can  be   imbedded   in  a   larger 
system  which  is  easier  to  solve.      The   larger  problem  can  for  example  be  a  discrete 
elliptic  problem  on  a   region  for  which  a   fast  Poisson  solver  can  be  used.      These 
techniques  are  also  of  interest  if  a   series  of  problems,  e.g.    the  same  elliptic 
equation  on  different  regions,   can  be   imbedded   in  the  same   larger  problem.    In  such 
a  case,   possibly  extensive  preprocessing  of  the   larger  problem  might  pay  off   if 
a   sufficient  number  of  smaller  problems     are   to  be  solved.     There  are  a  number  of 
different  versions  of  these  so  called  capacitance  matrix  methods.     Here  we  will 
only  consider  iterative  variants;   see  Dryja    (9),   O'Leary     and  Widlund    [17],     and 
Proskurowski   and     Widlund    (20,21).     We  also   limit  our  discussion  to  self  adjoint 
second  order  elliptic  problems  defined  on  a   region  which  is   imbedded   in  a   larger 
region. 
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We  can  then  adopt  the  same  notations  as  in  section  2,  and  regard  r3  as  the  boundary 
of  ni ,  which  is  the  region  of  interest.  The  complement  of  Q-]  u  r-,  with  resoect 
to  tne  larger  region  is  denoted  by  fjo.  The  vectors  xp  X2  and  'X3  are  associ- 
ated with  the  nodal  values  in  n, ,  n-  ^nd  r,  respectively.  To  simplify  our 
arguments,  we  assume  that  all  of  the  problems  have 'Symmetric ,  positive  definite 
coefficient  matrices.  ''-'.■. 


For  a  Neumann  problem  on  rii,  we  thfn  obtain,- 


Ml   "U. 

^3   "^33 
while  the  larger  problem  has  the  form  • 
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(6.2) 


The  value  of  b^  is  irrelevant  and  lye  can-  also,,  at  the  expense  of  one  solution  of 
equation  (6.2),  reduce  the  right  hand  si^e  of  (6.1).  tQ^the  c^ne  when  b-|  =  0.  Since 
we  can  solve  (6.2)  in  each  iteration  =s'te(5,  we 'can 'u'^e'the  operator  S  as  a  precondi- 
tioner.  Given  63  =  y  and  with  bi  and.  frj  2ero,  We  xSn  coispute  the  solution  of  equa- 
tion (6.2)  and  substitut.?  xi-afld  x,  into  •«<}ua-ttpn  (6.H;?..Wa  v.hen  obtain  the 

(11-1  '  ^■>--'3>'^.^ 

residual   b,  -  S       S'  y.     We  know  from^  section  3  that'thfi  operator  has  good 
spectral   properties.  ►       ■  ■   ■         :  ■' 


A  Dirichlet  problem  on  U]   can  be  reduced  tp  a  Neuinann  problem  on  Q2  by  the   follow- 
ing device.     The  finite  element  Dirichlet  problem   has  the  form, 

'^11   *^13  ]   "l  ]  . 

.°'  i';i1'*3l'' 

The  same  solution  can  also  be  obtained  from,  . 


Ml 


22 

T 
23 


M3 
^23 
^33 


■^r 

X2 

= 

.  '^3  ' 

"1   ' 

23 

(2) 

33 

"3 

^3! 

(6.3) 


which  has  the  unique  solution  (x-j.O.b,)  .  Equation  (6.3)  can  now  be  solved  by  the 
same  technique  as  the  previous  problem  on  n,   since  it  essentially  is  a  Neumann 
problem  on  a^.      The  resulting  iteration  matrix  is  S^^'S"  . 

It  is  also  of  interest  to  allow  more  general  coefficient  matrices  in  equation  (6.2). 
These  issues  will  be  discussed  in  a  forthcoming  paper,  see  Widlund  [23]. 


REFERENCES 

[1]      Bjerstad,   P.    E.,   Numerical    Solution  of  the  Biharmonic   Equation,   Ph.D.   Thesis, 
Stanford  Univ.    (1980). 

[2]      Bjfirstad,   P.    E.,   SIAM  J.    Numerical   Analysis  20   (1983)   59-77. 

[3]      Bj^rstad,   P.    E.   ind  Widlund,   0.   B.,     Solving  elliptic  problems  on   regions 
partitioned   into  substructures,   in:   Birkhoff,   G.    and  Schoenstadt,   A.    (eds.). 
Elliptic  Problem  Solvers   (Academic  Press,   New  York,    1984). 

[4]      Bj^rstad,   P.    E.   and  Widlund,  0.   B.,   Extended 'version  of  reference  3, 
in  preparation. 

[51      Ciarlet,   P.   G.,     The  finite  element  method  for_elliptic  problems   (North 
Holland,  Amsterdam,   1978J.  ^,~  ^^,  ■ 

[6]      Concus,   P.,   Golub,  G. ,   and  O'Leary,   D.,     Proc.    Symp.    Sparse  Matrix  Computa- 
tions,  Argonne  National    Lab.    (Academic  Press,'  19?5-). 

[71      Cottle,  R.,  Lin.   Alg.   Appl.   8  (1974)   189-211-. 

[8]      Dryja,  M.,   Numer.   Math.    39   (1982)   51-69. 

[9]      Dryja,  M.,     A  finite  element  -  capacitance -matrix  jnethod  for  elliptic  problems, 
to  appear. 

'      •  <■  ''?■". 

[101   Glowinski,   R.,   PerTaux,   J.   and  Dihn,'  Q.V.,  Domain  decomposition  methods   for 

nonlinear  problem?"  in-fTuid  dynamic's,-  IW^IA  *epo1-t^  No.    147   (July  1982). 

[Ill   Glowinski,   R.-i&^^.pi^qafie|^,,5^.,"'^S/$iJl_ev!4fw^  167-212. 

[121   Golub,  G.   H.  andr-Mayecs.  ^.■,  Tha'^  Proce'ettlng^.''' ' 

[13]    Grisvard,   P.,   Bosadary 'vac>lue'pr(S)V*ms--in  nonsni60th'   tlbmains,   Dept.   of  Math., 
Univ.   of  Maryland,^  Le^ture^otes  ,1^..(19§J3J  .^   ^.^,,.    ^,. 

[141    Hestenes,  M.,   Proc.    Symp.    Appl.   Math.   VI    (McGraw-Hill,   New  York,    1956)  83-102. 

[15]    Lions,  J.    L.   and  Magenes,   E.,   Nonhomo^eneous   boundary  value  problems  and 
applications,    I      (SpriTiger';'';fftlw  York,' 1972)  . .      ,'^    \ 

[16]    Luenberger,   D.    G.,    Introduction  to  linear  and  nonlinear  programming, 
(Addison-Wesley,    1973).      ■  ,     ...  "     i    ,     < 

[17]    O'Leary,   D.   and  Widlund,   0.,  "Math.   €91116.   33   (1979)   849-879. 

[181    Prezemieniecki ,   J.    5.,   Am.    Inst.   ^ero.   Astro.    7   (1963)    138-147. 

[19]    Prezemieniecki ,   J.    S.,   Theory  of  matrix   structural   analysis    (McGraw-Hill, 
New  York,   1968). 

[20]    Proskurowski,  W.   and  Widlund,   0.,   Math.   Comp.    30  (1976)   433-468. 

[21]    Proskurowski,  W.   and  Widlund.  0,,   SIAM  J.    Sci.    Stat.   Comput.    1    (1980)   410-468. 

[22]    Strang,  G.   and  Fix,   G.    J.,  An  analysis  of  the  finite  element  method, 
(Prentice-Hall,   1973). 

[23]    Widlund,   0.,   On   the  convergence  of   finite  element-  capacitance  matrix  methods, 
in  preparation. 

ACKNOWLEDGEMENTS 

This  work  was   supported   by   the  National    Science   Foundation  under  Contract  NSF- 
MCS-8203236     and  by   the  U.    S.    Department  of   Energy  under   contract   DE-AC02-76- 
ER03077-V     at  the  Courant  Mathematics  and  Computing  Laboratory. 


13 


This  book  may  be  kept 

FOURTEEN    DAYS 


A  fine  will  be  charged  for  each  day  the  book  is  kept  overtime. 

1 

CAY LORD  U2 

pnrNTED  IN  U   S  A, 

NY J  C.2 

Comp.  Sci.  Dept . 
TR-101   Widlund 
Iterative  methods  for 
elliptic  problems  on  ... 


LIBRARY 

N.Y.U.  Courant  Institute  of 

Mathematical  Sciences 

251  Mercer  St. 
New  York,  N.  Y.    10012 


'rm^^m^^ 


